Strongly localised molecular orbitals for a-quartz 

Oleh Danyliv § and Lev Kantorovich 

Department of Physics, Kings College London, Strand, London WC2R 2LS, UK 

Abstract. A previously proposed computational procedure for constructing a set 
of nonorthogonal strongly localised one-electron molecular orbitals (O. Danyliv, L. 
Kantorovich - Phys. Rev. B, 2004, to be published) is applied to a perfect a- 
quartz crystal characterised by an intermediate type of chemical bonding. The orbitals 
are constructed by applying various localisation methods to canonical Hartree-Fock 
orbitals calculated for a succession of finite molecular clusters of increased size with 
appropriate boundary conditions. The calculated orbitals span the same occupied 
Fock space as the canonical HF solutions, but have an advantage of reflecting the true 
chemical nature of the bonding in the system. The applicability of several localisation 
techniques as well as of a number of possible choices of localisation regions (structure 
elements) are discussed for this system in detail. 



PACS numbers: 31.15.Ar, 71.15.Ap, 71.20.Nr 



Submitted to: J. Phys.: Condens. Matter 



§ On leave from Institute for Condensed Matter Physics, National Academy of Science of Ukraine, 
Ukraine (e-mail: oleh.danyliv@kcl.ac.uk) 



Strongly localised molecular orbitals 



2 



1. Introduction 

A quantum cluster embedding [TJ El El HI El IE] has become a powerful computational 
tool in electronic structure theory of extended systems, such as large biological molecules 
[ll El El El E] surface defects and adsorption on crystal surfaces [101 IHJ E2j or points 
defects in the bulk of crystalline [T3*] ITi] or amorphous ^H] systems. The embedding 
methods originate from a model in which a single local perturbation is considered in the 
direct space of the entire system inside of a finite quantum molecular cluster in great 
detail, whereas a more approximate method is used to account for the rest of the system 
surrounding the cluster. 

A rather general embedding method is presently being developed in our laboratory. 
The central idea of our method is based on the exact partitioning of the entire system 
electron density into two components, one localised within the cluster and the other - 
outside it, i.e in the environment. Construction of overlapping (not orthogonal) strongly 
localised molecular orbitals (LMO's) as building blocks of the entire system is crucial 
for this technique. The LMO's are designed to represent the true electronic density of 
reference systems (such as e.g. 3D ideal perfect crystals or 2D periodic crystal surfaces) 
and are constructed to have transparent chemical meaning, e.g. to represent ions in the 
case of ionic systems and covalent bonds in covalently bound systems. Although we are 
not yet concerned in this study with biological systems which do not possess periodicity, 
we note that most of the ideas of our method can also be applied to these systems as 
well. 

A convenient and simple method for calculating the LMO's for 3D periodic systems 
(e.g. perfect crystals) was recently suggested in [Tfi] . This method is based on finding 
appropriate linear combinations of the canonical Hartree-Fock (HF) solutions for a 
sequence of finite molecular clusters of increased size. The linear combinations are 
chosen in such a way as to optimise special localising junctionals constructed to obtain 
orbitals localised within certain regions (e.g. bonds, atoms, ions, etc.). There may be 
several different regions in every unit cell of the crystal. 

The method was successfully applied in [HI] to two cases of extreme ionic (MgO 
crystal) and covalent (Si crystal) bonding. In both cases four LMO's were found in 
every unit cell. In the former case every unit cell is composed of a single region which 
is associated with an oxygen ion; the region contains eight electrons and is described 
by four mutually orthogonal LMO's. In the latter case (the Si crystal) every unit cell 
is represented by four neighbouring regions. Each region is associated with a pair of 
nearest Si atoms, contains 2 electrons and is described by a single double occupied LMO. 
The four regions belonging to the same unit cell have one common Si atom at the centre 
of a tetrahedron and other four Si atoms form its vertices. The four LMO's within the 
same cell do overlap and thus are not orthogonal. 

The main purpose of this paper is to check if our method [THj can be extended to 
systems which have more complicated types of chemical bonding. This is invaluable for 
the future development of our embedding method towards describing insulating crystals 
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with arbitrary type of bonding. Therefore, we here consider in detail the a-quartz 
(SiC>2) crystal, which may be thought of as a prototype system with an intermediate 
(ionic-covalent) type of chemical bonding. Essentially two main questions are addressed 
here with respect to the localisation of the calculated LMO's: (i) the choice of regions 
and (ii) the choice of localisation methods (i.e. the localisation functionals). 

The plan of the paper is the following. In the next section the main ideas of 
our method are briefly described (for the full discussion, see [H^)- Three localisation 
methods are introduced (one of them was not used in our previous study jTHI) alongside 
with a choice of three localisation criteria. Application of our method to the a-quartz 
crystal is considered in Sectional Brief conclusions are given in Section 0] 

2. Localisation procedure 

2.1. General approach 

Let us assume that we know an occupied canonical set of one-electron molecular orbitals, 
{<fi(r)}, for a perfect 3D periodic crystal. These orbitals may be obtained as eigenvectors 
of the appropriate Hartree-Fock (HF) problem using e.g. the CRYSTAL code which 
employs directly the periodic symmetry. In our method, however, we consider instead 
a specially designed set of finite clusters of increased size and find the HF solutions 
for them using one of the available quantum chemistry packages. It was demonstrated 
in [TE] that this approach is equivalent to using a periodic-crystal electronic structure 
approach as far as the LMO's are concerned, provided that large enough molecular 
clusters are used. 

The canonical molecular orbitals (CMO's) are orthogonal and span the entire 
occupied Fock space. They are not localised in space, and have a non-zero contribution 
on atoms of every unit cell in the crystal. In practice, when the cluster method is 
employed, they span the entire cluster. In other words, the CMO's are assumed to be 
given as a linear combination of the atomic orbitals, X;u( r )> centred on all atoms of the 
cluster in question: 



In order to describe the crystal as a set of overlapping (non-orthogonal) localised 
functions, {(p a (r)}, which span the same occupied Fock space, one has to obtain 
appropriate linear combinations of the original canonical set {^(r)}. In order to do 
this, it is first necessary to identify regions of space where each of the functions <£> Q (r) 
has to be localised. Although any (non-singular) linear combination of the canonical 
set will give the same electron density p(r), we adopt here a strategy based on the 
chemistry of the system in question. Namely, the choice of the localisation regions 
in the first instance is based on the expected type of the chemical bonding in the 
system, e.g. atoms/ions in the cases of atomic/ionic systems, two nearest atoms in 
the case of covalent bonding, etc. A more complicated choice is anticipated in the 




(1) 
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cases of intermediate bonding as will be demonstrated in Section OJ Several different 
nonequivalent regions may be necessary to represent a crystal unit cell which can then 
be periodically translated to reproduce the whole infinite crystal. Note that several 
localised orbitals may be associated with each region. For instance, in the case of the 
Si crystal one needs four localised regions associated with four bonds; each bond is 
represented by a single double occupied localised orbital and all four bonds have one 
common Si atom in the centre of the tetrahedron. 

Once the localised regions are identified, it is necessary to find linear combination 
of the CMO's which are localised in each of the regions, 

occ 

&(r) = £ U aj ^(r) = £ ^(r) (2) 

j A* 

The transformation U =||C/ a j|| of the CMO's within the occupied subspace is arbitrary 
and, in general, non-unitary. In the latter case the expression for the density via the 
new set of orbitals should contain the inverse of the overlap matrix S = ||S a &|| [TH] : 

occ ^ _ ^ 

p(r) = 2££ a (r)(s) a6 ^(r) (3) 

ab 

where S a b = ((p a (r)\ ^fe( r )) is the overlap integral. The double summation here is 
performed over all localised orbitals of the whole infinite crystal. If the transformation is 
unitary, then both the overlap matrix and its inverse are unity matrices and the density 
takes on its usual "diagonal" form. 

In general, any localisation procedure is equivalent to some transformation U of 
the CMO's. To find the necessary transformation for, say, region A, an optimisation 
(minimisation or maximisation) problem is formulated for some specific localising 
functional Qa [{<Pa}] with the constraint that the LMO's associated with region A are 
orthonormal. This leads to a standard eigenvalue-eigenvector problem: 

occ 

J2^paj = KUai (4) 
3 

for the elements of the transformation matrix U. Here is a matrix element of an 
operator VLa calculated using canonical orbitals v?i( r )- The operator Qa is uniquely 
defined from the functional [{<£>„}] . Although for some localising functionals (see, 
e.g. ^Hj) the matrix elements Qfj may depend on the LMO's themselves so that the 
problem (jlj) is to be solved self-consistently, we do not consider those functionals in this 
paper. 

Note, that LMO's associated with different regions will not be orthogonal in this 
method. This is because they are obtained from different localising functionals which 
strongly depend on the region in question, so that LMO's from different regions are 
determined by solving different secular problems. For instance, if LMO's {(p a (r)} 
correspond to region A, then the LMO's {^(r) = y5 (r — L)} are obtained for a 
physically equivalent region B separated from A by a lattice translation L. 

Using a physical argument, each region is associated with a certain even number 
of electrons 2n. Therefore, if VLa is minimised, we choose the first n eigenvectors of the 



Strongly localised molecular orbitals 



5 



problem (J3J); if, however, Qa is maximised, the last n solutions are adopted. By collecting 
LMO's from all regions in the unit cell and then translating those over the whole crystal 
it should be possible to span the whole occupied Fock space and thus construct the total 
electron density 0- The larger the finite cluster used while calculating the canonical 
orbitals, the closer the Fock space will be reproduced by the LMO's. 

To summarise, we first suggest a possible set of localisation regions in the unit 
cell and then consider a set of finite molecular clusters (with appropriate boundary 
conditions) which have all these regions in their central part. Then, we obtain the 
occupied canonical HF orbitals for each of the clusters using a standard quantum- 
chemistry package. Out of all the clusters considered, a cluster is chosen for which 
the electron density is well converged in its central part. Next, using a localisation 
functional, canonical occupied HF orbitals of the chosen cluster are transformed into 
LMO's. The procedure is repeated for several localisation functionals and in each case 
localisation criteria are applied. Then, if necessary, a different choice of localisation 
regions is made, and the whole procedure is repeated. As will be seen in Section El in 
the case of the Si0 2 crystal, three different sets of regions can be suggested; however, 
the same set of clusters will be used to calculate the LMO's in each case. 



2.2. Localising functionals 

In a number of methods ^J] the localising functionals are proportional to the non- 
diagonal electron "density" associated with region A, 

n 

(t a (t,t') = £<£ a (r)#;(r / ) (5) 

a=l 

where the summation is performed over all n LMO's of region A. Note that for 
convenience of the final equations we have omitted a factor of two above as it is 
unimportant for the eigenproblem (j3J) to be solved. Therefore, the functionals can 
be represented in the following general form: 

Ti Tl OCC 

V A = ! [*Wr,rO] . * = £ / ft(r)*W.(r)dr = £ £ U^n^U*® 

" ^ — 1 ** * — 1 ,-7„ 



a=l a=l jk 



where ft a is the localisation operator and the Hermitian matrix Q A = \\ Q A k || can easily 
be written in terms of the canonical MO's using the definition (fT|): 

Qf k = U q a Wl) = £ c%c c vk ( X ,\ n A |x„> (7) 



For all methods to be considered below both the operator Qa and the matrix £l A do not 
depend on the LMO's sought for, so that in order to obtain the localised orbitals one 
has simply to find the eigenvectors of the matrix fl A using Eq. (J3J). Three particular 
localisation methods implemented in this work are considered in the following in more 
detail. Note that one of the methods (method G) was not considered in p| 
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2.2.1. Mulliken's net population (method M) In this method the localisation region A 
is specified by a selection of AO's (e.g. on one or two particular atoms in the unit cell). 
Then, the net atomic Mulliken's [201 population produced by the LMO's in the selected 
region is maximised [2B EH E2 • in this case 

^fk = E C^jS^ w C^ k (8) 

where is the overlap integral between two AO's x M and Xu- The summation here 
is performed over AO's which are centred in the chosen region A. This way one can 
make the LMO's to have a maximum contribution from the specified AO's in region A. 
Sometimes a different choice of AO's centred on the same atoms may lead to physically 
identical localisation; however, this is not the case in general [Ej. This method will be 
referred to as method M. 

2.2.2. Mulliken's gross population (method G) If, instead, the Mulliken's gross 
population on the atoms belonging to region A is maximised, one arrives into the Pipek- 
Mezey localisation scheme [221 OH] ■ I n this case the expression for fl^ k is very similar to 
that given by Eq. (JBJ): 

^fk = 9 E E \p^^vk + C^jSvpC^A (9) 

The first summation here is performed over AO's which are centred in the chosen region 
A and another summation is performed over all AO's. This method will be referred to 
as method G. 

2.2.3. The projection on the atomic subspace (method P) The Roby's population 
maximisation j2H| gives LMO's for which the projection on the subspace spanned by 
the basis orbitals centred within the selected region A is a maximum, or is at least 
stationary PTJJ EI] . In this method the localisation operator fl^ in Eq. (jUJ) is chosen in 
the form of a projection operator, so that: 

(10) 

where S^ 1 stands for the inverse of the overlap matrix defined on all AO's /i, v G A. 
Here the first double summation is performed over all AO's of the system. Note that 
the idempotent operator Qa projects any orbital into a subspace spanned by the AO's 
associated with region A only. Therefore, by choosing particular AO's (and thus the 
region) one ensures the maximum overlap of the LMO's with them. It is seen that this 
method, which will be referred to as method P, although different in implementation, is 
very similar in spirit to the previous two methods. 



E 
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2.3. Localisation criteria 

An application of the various schemes described above results in LMO's which are 
localised in the 3D space differently. It is therefore useful to have simple criteria which 
can identify the degree of their localisation. Note that each of the localisation methods of 
the previous subsection corresponds to a particular linear combination of the canonical 
orbitals and thus will result in exactly the same electron density (J3J) provided, of course, 
that a sufficiently large cluster has been used in the LMO's calculation. We assume in 
this section that this is always the case. 

Three methods will be used to assess the localisation of calculated LMO's jTH]. 



2.3.1. Localisation index The first method was proposed by Pipek and Mezey 
is based on the calculation of the so-called localisation index: 

-i 



d„ 



E 

B 



fi&B v 



and 



where the first summation is run over all atoms B of the entire system. Here the quantity 
in the square brackets is similar to the diagonal part of the localisation operator matrix 
(JHJ) calculated on real localised orbitals. Qualitatively, the localisation index gives the 
number of atoms on which the orbital <£> a (r) is predominantly localised. Therefore for 
the ionic type of bonding one would expect d a ~ 1 and for a valence LMO describing a 
covalent bonding d a ~ 2. 



2.3.2. Eigenvalues of the overlap matrix Alternatively, the overlap between LMO's 
also gives an important information about their localisation. That is why as the second 
criterion we shall consider the maximum eigenvalue of the overlap matrix. Note that 
for periodic structures it is more convenient to use the Fourier transformation of the 
overlap matrix [2*lj : 

S ab (k) = £(^(r)|^(r-L))e* kL (12) 

L 

where k is a point in the Brillouin zone, y5 a ( r ) and £>b(r) are LMO's in the zero (central) 
elementary unit cell and L is the lattice translation vector. Note, that if any of the 
eigenvalues, A(k), of the overlap matrix S = ||S' a 6(k)|| is found to be larger than 2, it is 
impossible to obtain the total crystal density in this basis by expanding the inverse of 
the overlap matrix in Eq. (JHJ) in powers of the overlap (the Lowdin's method, see j21] 
for a detailed discussion). Therefore, the existence of large eigenvalues of the matrix S 
corresponds to a weak localisation of the LMO's. 



2.3.3. Gap in eigenvalues of the localisation problem The eigenvalues A a of the secular 
problem (jlj) can also be used to indicate the degree of localisation [THj. Indeed, if the 
localisation functional Qa used is appropriate, then (i) the chosen n solutions would 
have close eigenvalues A a which correspond to their similar localisation in region A, 
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and (ii) the gap AA in the eigenvalues A a between the chosen n and other solutions 
is considerable, i.e. the other solutions will correspond to much worse localisation in 
region A (cf. (23)- Therefore, in order to check the localisation of the LMO's, we shall 
also use the parameter AA. 

One may assume that better localisation will give larger gap value AA. In principle, 
if the given region and the right number of LMO's n, associated with it, were chosen 
correctly, one should expect some gap AA in the eigenvalues A a . 



3. Si0 2 bulk 

a-quartz (Si02) crystal has a hexagonal Bravais lattice and corresponds to the _D| 
(P3i21, No. 152) space group symmetry [2B1I27]. The equilibrium crystal structure 
has been found using the Density Functional Theory, plane wave basis set, periodic 
boundary conditions and the method of ultrasoft pseudopotentials as implemented in 
the VASP code j2HH23E0I- The calculation was carefully converged with respect to the k 
point sampling and the plane wave cut-off. The lattice found is specified by elementary 
translations a± = a(0, — 1,0), a 2 = a(|,^,0) and a3 = c(0,0, 1) with a=4.913 A 
and c=5.4046 A. Each cell contains nine atoms distributed over three Si0 2 molecules. 
Three 3a positions (the local symmetry C 2 ) are occupied by Si atoms, the position of 
the first Si atom is given by the fractional coordinates (—u, —u, |) with w=0.4697; six 
oxygen atoms occupy a general position 6c which can be generated from the fractional 
coordinates (x,y,z) using x=0.4135, y=0.2669 and 2=0.1191. The tetrahedron of 
oxygen atoms with a silicon atom in its centre is almost regular with slightly different 
Si-0 distances of J2(Sil-0)=1.613 A and fl(Si2-O)=1.604 A. The whole 3D crystal 
structure can be constructed from SiiOSii units connected together at the Si atoms as 

4 4 

shown schematically in Fig. [T^a). Four such units have a common Si atom at the centre 
of a tetrahedron, six complete inequivalent units (positioned differently in space) form 
an elementary cell. Note that this is similar to the Si crystal structure where the whole 
crystal can be composed of Si i Si i units [To] . 

Because of such an arrangement of atoms in the a-quartz crystal, it is reasonable 
to assume that every unit SiiOSii forms one independent region. However, as will 

4 4 

be shown in the following, one can alternatively consider two or three regions made 
out of each unit as well. Therefore, in choosing a set of finite clusters, we ensured 
that the whole SiiOSii unit was in the centre of each of the quantum clusters. Three 

4 4 

clusters were considered: Si 2 7 , Si 8 25 and Si 40 Oio3, containing 9, 33 and 143 atoms, 
respectively. The smallest and middle size clusters used are shown in Fig. ^ To create 
proper termination of the clusters at their boundary, we implemented special boundary 
conditions suggested in [THj. In particular, we surrounded clusters by pseudoatoms Si* 
which are directly connected to the cluster O atoms. Each Si* pseudoatom is made of 
"classical (3/4-th) and "quantum" (1/4-th) parts. The "classical" part is represented by 
the electrostatic potential due to a +1.8e point charge (which is a 3/4-th of the effective 
charge on a Si atom in the lattice), e being the elementary charge. The "quantum" part 
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of the Si* pseudoatom consists of a central repulsive electronic potential V(r) (added to 
mimic the screening of the Si core by the valence electrons) and a single valence electron. 
The parameters of the potential and the basis set for the pseudoatoms were optimised 
to get proper effective charges on Si and O atoms and to eliminate the contribution 
of the Si* electron at the top of the valence and the bottom of the conduction bands. 
Note that the described boundary conditions were found to be crucial only for the 
smallest cluster; for other two clusters simpler boundary conditions (e.g. termination 
by hydrogen atoms) were also tried and found to give practically identical results for the 
LMO's and thus will not be discussed further. To simulate the Madelung field, clusters 
were surrounded by an array of nearly 2.5 x 10 4 point charges, containing +2.4e charges 
to mimic Si atoms and -1.2e charges for oxygens. 

To simplify the initial HF calculations required to check the convergence of the 
electron density with the cluster size and generate all occupied canonical molecular 
orbitals, only valence (3s 2 3p 2 on Si and 2s 2 2p 4 on O) electrons were considered explicitly 
by using for both species the coreless HF pseudopotentials (CHF) with LP-31G basis 
set from Ref . [3T] ■ The calculations were made with the use of the Gamess-UK package 

021 . 

Fig. El shows the convergence of the HF electron density with the size of the cluster 
along the Sil-0 direction. We carefully checked, by plotting the densities along other 
directions and by making 2D plots, that this direction is representative for assessing 
the convergence in this system. The deep minimum in the density at the oxygen atom 
is due to the pseudopotential method used. One can see that the difference between 
curves for the middle sized and the largest clusters is negligible. This suggests that the 
middle sized cluster SigO^ is perfectly sufficient for our purposes and thus was used in 
all the calculations described below. 

As has been mentioned above, the elementary "brick" we can build the system from 
is the unit SiiOSii shown in the centre of each of the clusters in FkrH] 

cLS cL Si 1-0-S12 

4 4 

molecule. Each such unit should be assigned 8 electrons in total: 6 electrons come from 
the O atom and by 1 electron from each of the two Si atoms. Note that each Si atom 
contributes to four different units which contain this Si atom. This simple analysis 
allows us to suggest at least three possible choices (models) for the localisation regions: 

(i) Each unit SiiOSii containing 8 electrons is considered as a single region; therefore, 

4 4 

to obtain the corresponding four (double occupied) LMO's, one has to choose all 
AO's centred on the atoms Sil, Si2 and O in the centre of the cluster when applying 
any of the localising functionals discussed above. Note that all four LMO's obtained 
using this partition method will be orthonormal as eigenvectors of the same secular 
problem (j3J). 

(ii) Each pair of atoms Sil-0 and Si2-0 can be considered as a separate region, i.e. 
there will be two regions in total to describe every unit SiiOSii; four electrons 

4 4 

distributed over two (double occupied) LMO's will be associated in this case with 
each of the two regions. Two LMO's associated with either of the two regions will 
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Figure 1. The first two quantum clusters, Si2C>7 (a) and SisC^s (b), used in our 
HF calculations. Point charges surrounding the clusters are not shown. Both silicon 
atoms (Sil and Si2) and the oxygen atom of the central unit SiiOSu are indicated in 
each case. The elementary "brick" the whole system can be composed from is shown 
schematically in (a). 

be mutually orthogonal; however, the LMO's belonging to different regions will 
have a non-zero overlap. 

(iii) Finally, each unit can be split into three different regions: (i) the first region, 
containing two electrons and described by a single LMO, is constructed to describe 
a covalent bond Sil-O; this can be achieved by enforcing localisation on 2p AO's of 
the O atom and all AO's of the Sil atom; (ii) the second region is formed similarly 
to describe the Si2-0 bond; (iii) finally, the remaining four electrons are attached 
to the third region which is localised predominantly on the O atom giving rise to 
two more (double occupied) LMO's; in this case the 2s AO's centred on the O atom 
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1 .2 




distance (A) 

Figure 2. The HF electron densities for the three clusters are plotted along the Sil-0 
direction. The convergence is obvious: the curves for the middle sized and the largest 
clusters are indistinguishable. 

can be used to inforce localisation. Thus, in this case there will be three sets of the 
LMO's: two orthogonal LMO's belonging to the O region and other two LMO's 
belonging to the "bond" regions; the latter two LMO's have a non-zero overlap 
with any other LMO. 

For each choice of the localisation regions described above (which will be referred to 
as regions models hereafter), we can apply either of the three localisation methods of 
Section El to calculate the LMO's using the obtained occupied canonical orbitals of 
the middle cluster. When solving the secular problem the contributions of the 
boundary pseudoatoms Si* were removed from the canonical orbitals which then were 
renormalised. Using the obtained LMO's, one should appropriately translate and rotate 
them in order to obtain all LMO's comprising the whole primitive cell. (For instance, 
there will be six sets with four LMO's in each in the primitive cell for model 1.) By 
applying lattice translations to all the LMO's associated with the primitive cell, the 
whole infinite crystal is reproduced. It is then possible to calculate the total electron 
density p(r) of the whole crystal using Eq. (JH1). The necessary lattice summations are 
handled exactly by converting into the k space [21] . These calculations have been done 
for all nine cases (three localisation methods versus three choices of the localisation 
regions). The calculated p(r) matched exactly the original HF density in the central 
part of the cluster in all cases indicating that a very good degree of localisation was 
achieved in each case. As an example, a 2D contour plot of the total valence electron 
density in the plane of the molecule Sil-0-Si2 calculated using LMO's obtained by 
method M in model 1 is shown in Fig. 01 The cross-section of this plot corresponds to 
a solid line (the middle cluster) in Fig. El One can see that a considerable amount of 
the charge is concentrated around oxygen atoms. 




Figure 3. The total electron valence density p(r) in the Sil-0-Si2 plane by combining 
the contributions of LMO's (method M, model 1) across the infinite crystal. The 
density was calculated using the method TIj based on Eq. ©. 



The partial region densities, pa(t) = 0^(r, r) (see Eq. (JHJ)), corresponding to each 
of the regions and generated from the LMO's calculated using method M, are shown in 
Fig. 0] as closed 3D surfaces of constant density. The value of the density is chosen in 
such a way so that 90% of the region electron charge be contained inside every surface. 
In the case of model 1 (the upper panel) only one density is shown; in the case of model 
2 (the middle panel) two densities are shown simultaneously, while in the third case 
(model 3, the bottom panel), all three partial densities are presented. One can see that 
in the cases of models 2 and 3 the LMO's belonging to neighbouring regions strongly 
overlap. As was mentioned before, the LMO's belonging to different regions for these 
two models are not orthogonal. At the same time, the overall shapes of the density for 
each of the region models are very similar demonstrating a clear aggregation around the 
O atom in the middle of the SiiOSii unit in agreement with the total density shown in 

Fig. m 

A comparison of the LMO's calculated using three localisation methods is presented 
in Fig. El for region model 1. In this figure, the partial densities pa{t) are shown in 
each case along the Sil-0 direction through the SiiOSii unit. It is clear that the 

4 4 

partial density obtained from LMO's calculated using method M are found slightly 
more localised, whereas the localisation obtained using method P is slightly worse than 
given by the two others. Still, the difference between the densities is extremely small so 
that we can conclude at this point that all three techniques perform practically equally 
well (at least for the system under discussion). 

The picture becomes more complicated, however, at least at the first sight, when 
the localisation criteria of Section 12.31 are applied in each of the nine cases as shown 
in Table ^ Three rows in this table correspond to the three different region models; 




Figure 4. Constant partial density plots for three choices of localisation regions: (a) 
model 1, when the whole unit Sii OSij. is associated with a single region; (b) model 2, 

4 4 

when two regions Sil-0 and Si2-0 are identified and (c) model 3, when three regions, 
Sil-O, Si2-0 and O, are identified. In every case the value of the density shown is 
chosen to enclose 90% of the total electron charge associated with the region. The 
localisation method M was used in each case. 



in each case there are four LMO's in total which are occupied by eight electrons of 
the elementary SiiOSii unit. The three columns in the Table correspond to the three 
different localisation methods used, and each of the localisation criteria is shown for 
every method. 

The fist criterion (the localisation index d a ) is slightly above one in most cases, and 
is smaller than two in all nine cases. This means that the LMO's are mostly localised 
on a single atom with some contribution coming from the nearest atoms. However, 
the maximum eigenvalues of the overlap matrix were found to be below two only for 
the region model 1 and localisation methods M and G; in the cases of models 2 and 
3 eigenvalues around three were found indicating worse localisation. Finally, the gap 
between the eigenvalues of the secular problem of Eq. (J1J, A A, was found to be too 
small for the region models 2 and 3, whereas in the case of model 1 and for all three 
localisation methods the gap is considerable, especially for methods M and G. This 
means that the choice of the regions in models 2 and 3 is somewhat artificial which is 
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distance (A) 



Figure 5. The partial electron density for region model 1 calculated using different 
localisation methods and plotted along the Sil-0 direction. 



not surprising because of a very strong overlap between LMO's corresponding to the 
neighbouring regions in these two models. 

It follows from this analysis that the model 1, in which the whole elementary 
unit SiiOSii is considered as one region, results in the best localisation of the LMO's, 

4 4 

especially if methods M and G are used. 

In spite of subtle differences in the applied localisation criteria which seem to favour 
the model 1 and the methods M and G, we stress that very good localisation of the 
LMO's was obtained in all cases. This conclusion is also supported by an observation 
that LMO's generated within different models (choices of the regions) span the same 
Fock space. To make such a conclusion, we calculated the projection of the LMO's of 
models 2 and 3 on the space spanned by the four orthogonal LMO's obtained in model 
1 and then subtracted the projection from the original orbitals. The calculated residual 
parts were found negligible in all cases. Note, that this conclusion is not obvious because 
each LMO depends on all AO's of the entire cluster. 

4. Conclusions 

In this paper we have calculated strongly localised molecular orbitals (LMO's) for the 
Si0 2 crystal (a-quartz) using the method developed earlier [TH]. The starting point for 
the choice of the localisation regions was an observation that the whole crystal can be 
reproduced by rotating and translating a single elementary unit SiiOSii , containing an 

4 4 

O atom and quarters of the two Si atoms which the O atom is directly connected to. 
Three localisation methods were applied and three models for choosing the localisation 
regions were tried in each case: (1) the whole unit was considered as one region; (2) 
the unit was split into two and (3) three regions. Although in all cases well localised 
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Table 1. Localisation criteria calculated for all region models and localisation 
methods. 

orbitals were obtained, we find that the first choice of the localisation region in which 
the whole unit was chosen as a single region, is preferable. 

The LMO's produced in models 2 and 3 were found very close to a linear 
combination of the orthonormal LMO's obtained within model 1. If taken from 
the nearest regions belonging to the same elementary unit, they appeared to have a 
significant overlap with each other. On the other hand, the LMO's belonging to different 
units (in either of the models) do not overlap strongly which is confirmed by various 
localisation criteria applied in this work and by the corresponding plots of the partial 
densities. 

Since our previous calculations reported in ^H] were done for the extreme cases of 
ionic (MgO) and covalent (Si) bonding, it follows from the results of the present study 
that our method is also applicable to the crystals with intermediate types of chemical 
bonding. 

Note that we did not consider in this study localisation method E ^H] based on 
the energy minimisation of the structure element corresponding to the chosen region. 
This is because it was found in ^B] that the orbitals obtained by this method for the Si 
crystal were not sufficiently localised. 

Although the LMO's reported in this work may be useful to characterise the 
chemical bonding in the given crystal, they are needed for the embedding method which 
is under development. Different possibilities in choosing localisation regions open up 
various ways in terminating the quantum cluster when considering, e.g. a point defect 
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in the crystal bulk or an adsorbed species on the crystal surface. This variety of options 
may be extremely useful in applications to keep the size of the cluster as small as 
possible. If, for instance, one would like to terminate the cluster with Si atoms, then 
either of the region models can be used (model 1 would probably be more convenient 
as the orbitals within each region are orthonormal). However, if a termination with 
oxygens is required, then region models 2 or 3 may prove to be more useful. In practice, 
a combination of terminations may be preferable, when both Si and O atoms are used 
at the boundary. In those cases all three models for choosing localisation regions may 
be employed. 
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